SECOND-ORDER HYPERBOLIC FUCHSIAN SYSTEMS. 
GOWDY SPACETIMES AND THE FUCHSIAN NUMERICAL ALGORITHM 
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Abstract. This is the second part of a series devoted to the singular initial value problem for 
second-order hyperbolic Fuchsian systems. In the first part, we defined and investigated this 
general class of systems, and we established a well-posedness theory in weighted Sobolev spaces. 
This theory is applied here to the vacuum Einstein equations for Gowdy spacetimes admitting, by 
definition, two Killing fields satisfying certain geometric conditions. We recover, by more direct 
and simpler arguments, the well-posedness results established earlier by Rendall and collaborators. 
In addition, in this paper we introduce a natural approximation scheme, which we refer to as the 
Fuchsian numerical algorithm and is directly motivated by our general theory. This algorithm 
provides highly accurate, numerical approximations of the solution to the singular initial value 
problem. In particular, for the class of Gowdy spacetimes under consideration, various numerical 
experiments are presented which show the interest and efficiency of the proposed method. Finally, 
as an application, we numerically construct Gowdy spacetimes containing a smooth, incomplete, 
non-compact Cauchy horizon. 

1. Introduction 

This is the second part of a series [31 H] devoted to the initial value problem associated with 
the Einstein equations for spacetimes endowed with symmetries. We are typically interested in 
spacetimes with Gowdy symmetry, and in formulating the Einstein equations with data on a singular 
hypersurface, where curvature generically blows-up. In the first part [3 , we introduced a class of 
partial differential equations, referred as as second-order Fuchsian systems, and we established a 
general well-posedness theory within Sobolev spaces with weight (on the coordinate singularity). 
In the present paper, we tackle the treatment of actual models derived from the Einstein equations 
when suitable symmetry assumptions and gauge choices are made. 

We consider here (3-1- l)-dimensional, vacuum spacetimes {M,g) with spacelike slices diffeomor- 
phic to the torus T^, satisfying the vacuum Einstein equations and the Gowdy symmetry assump- 
tion. That is, we assume the existence of an Abelian T^-isometry group with spacelike orbits and 
vanishing twist constants [S] . These so-called Gowdy spacetimes on were first studied in . In 
the past years, a combination of theoretical and numerical works led to a detailed understanding 
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of the Gowdy spacetimes, achieved by analyzing solutions to the Einstein equations as the singular 
boundary is approached; cf. [13 UHl [201 [HI [22 ■ 

Before we can indicate precisely our contribution in the present paper, let us provide some 
background on Gowdy spacetimes. Introduce coordinates {t, x, y, z) such that (x, y, z) describe 
spatial sections diffeomorphic to while i is a timelike variable. We can arrange that the Killing 
fields associated with the Gowdy symmetry coincides with the coordinate vector fields in a 

global manner, so that the spacetime metric reads 



(1.1) 



3=4= e^^^{-dt^ + dx^) + i {e^{dy + Qdzf + e^^dz), t > 0. 



Hence, the metric depends upon three coefficients P = P{t,x), Q = Q{t,x), and A = A{t,x). We 
also assume spatial periodicity with periodicity domain U := [0,27r). 

In the chosen gauge, the Einstein's vacuum equations imply the following second-order wave 
equations for P, Q, 



(1.2) 



Ptt + f-P.. = e''^{Qt-Qi), 



Qt 



t 



Qxx 



-2{PtQt - P^Q^) 



which are decoupled from the wave equation satisfied by the third coefficient A: 
(1.3) Au - A,, - - P2 + e2P(g2 _ q2^^ 

Moreover, the Einstein equations also imply constraint equations, which read 
(1.4a) A, = 2t {P^Pt + e^^'Q^Qt) , 

(1.4b) At = t + te^PQl + p2 ^ g2 



It turns out that (|1.3|) can be ignored in the following sense. Given a time to > 0, we can 



prescribe initial data (P,Q)\tg for the system (1.2) while assuming the condition 



(1.5) 



„2P 



QxQt) dx = 



Sit t = tn 



Then, the first constraint in (1.4) determines the function A at the initial time, up to a constant 



which we henceforth fix. Next, one easily checks that the solution (P, Q) of (1.2) corresponding 
to these initial data does satisfy the compatibility condition associated with (1.4) and, hence, 
(1.4) determines A uniquely for all times of the evolution. Moreover, one checks that (1.3) is 



satisfied identically by the constructed solution (P,Q,A). Consequently, equations (1.2) represent 



the essential set of Einstein's field equations for Gowdy spacetimes. We refer to ( 1.2 1 as the Gowdy 



equations and focus our attention on them. One could also consider the alternative viewpoint which 
follows from the natural 3 + 1-splitting and treats the three equations (1.2)-(1.3| as an evolution 



system for the unknowns (P, Q, A), and ( 1.4 ) will be regarded as constraints that propagate if they 
hold on an initial hypersurface. 

Rendall and collaborators [H [THl [Hj developed the so-called Fuchsian method to handle singular 
evolution equations derived from the Einstein equations. This method allows one to derive precise 
information about the behavior of solutions near the singularity, which was a key step in the general 
proof of Penrose's strong cosmic conjecture eventually established by Ringstrom [22]. For the most 
recent developments we refer to 13 [Gj [7] ; especially, in [6] , a generalization of the standard Fuchsian 
theory was recently introduced. 
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In the first part [5], we revisited the Fuchsian theory and developed another approach to the 
well-posedness theory in Sobolev spaces which apphes to the singular initial value problem of a 
class of (second-order, hyperbolic) Fuchsian systems. The theory in [3] is particularly well-adapted 
to handle the Gowdy equations, as we show here. In fact, we propose here a rather simple proof of 
the well-posedness of the singular initial value problem for the Gowdy equations, which provides an 
alternative to the approach in Rendall [TH]. Moreover, in passing, we are also able to clarify certain 
issues. 

A second objective in the present paper is to present a new numerical approach and apply 
it specifically to Gowdy spacetimes. This approach is inspired by a pioneering work by Amorim, 
Bernardi, and LeFloch [T] , which computed solutions to the Gowdy equations with data imposed on 
the singularity. The approximation scheme proposed in the present paper, which we refer to as the 
Fuchsian numerical algorithm, is directly derived from our existence theory [3] and relies strongly 
on the hyperbolicity of the equations, so that error estimates for the numerical approximations are 
expected to hold. 

This paper is organized as follows. Section[2]is devoted to the theoretical discussion of the singular 
initial value problem for the Gowdy equations. First, we recall some heuristics for the behavior of 
the solutions at the "singularity" t = 0. Then, we show that our notion of singular initial value 
problems in [S] is consistent with these heuristics, and then state the well-posedness results that 
follows from the first part. In Section |4| we introduce a general numerical approach to solve second- 
order hyperbolic Fuchsian systems numerically and then apply it to the Gowdy equations. Several 
numerical experiments are presented and, finally we numerically construct a Gowdy symmetric 
solution of Einstein's field equations with a smooth non-compact Cauchy horizon. 



2. Singular initial value problem for the Gowdy equations 

2.1. Heuristics about the Gowdy equations. We provide here a formal discussion which mo- 
tivates the following (rigorous) analysis. Based on extensive numerical experiments, it was first 
conjectured (and later established rigorously) that as one approaches the singularity the spatial 
derivative of solutions {P,Q) to (1.2) become negligible and {P,Q), should approach a solution of 
the ordinary differential equations 

Q« + ^ = -2P.Q.. 

These equations are referred to, in the literature, as the velocity term dominated (VTD) equations. 
Interestingly enough, they admit a large class of solutions given explicitly by 

P{t,x)^\n{at''{l + Ct-^'')), 
(2.2) ^i-2fc 

where x plays simply the role of a parameter and a, C., ^, k are arbitrary 27r-periodic functions of x. 

Based on (2.2), it is a simple matter to determine the first two terms in the expansion of the 
function P near t — 0, that is, for fc 7^ at least, 

P(t x) 

lim^^^\imtPt{t,x)^-\k\, 
t-i-o Int t-i-o 
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hence 



lim(P(t,a:) + |fc(.T)|lni) = if{x), 



Similarly, for the function Q we obtain 



lim Q{t, x) = q{x), q := < ^ 



In a, k < 0, 

In(aC^), fc > 0. 



k<0, 
k^O, 

k>0, 



C 



A: < 0, 



0, fc = 0, 
^ fc > 0. 



(2.3) 



From (2.2), we thus have the expansion 

P = -|fc|lnt + ^ + o(l), 
Q = g + t2|^IV. + o(t2|fe|), 



in which fc, 1^9,(7, are functions of x. In general, P blows-up to +00 when one approaches the 
singularity, while Q remains bounded. Observe that the sign of fc is irrelevant as far the asymptotic 
expansion is concerned, and we are allowed to restrict attention to fc > 0. 



By plugging the explicit solution in the nonlinear terms arising in (2.1 1 one sees that e Qf is 



at 



of order i^d*^!^!) which is negligible since the left-hand side of the P-equation is of order t 
least when fc ^ 0. On the other hand, the nonlinear term PtQt is of order t^dfel-i)^ which is the 
same order as the left-hand side of the Q-equation. It is not negligible, but we observe that PtQt 
has the same behavior as —{\k\/t) Qt. 



In fact, observe that the homogeneous system deduced from (2.1) 



(2.4) 



t 



= 0, Qtt 



1 - 2fc 
t 



Qt -0, 



is solved precisely by the leading-order terms in (2.3 1. This tells us that, as t 
e^^Ql is negligible in the first equation in (2.1 ), while PtQt + (|fc|/i)<5t is negligible at t 



0, 



the term 
= 0. This 

discussion hence allows us to conclude that as far as the behavior at the coordinate singularity t = Q 



is concerned, the nonlinear VTD equations (2.1) are well approximated by the system (2.4) 



We return now to the nonlinear terms which were not included in the VTD equations, but yet 



are present in the full model (1.2). Allowing ourselves to differentiate the expansion (2.3), we get 



the following leading-order terms aX t — Q: 



„2P , 



2e^'^\kyj.ip\nt ^ 



qx = 0, 
= 0, 



|fcU = o, VxT^o, 
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' -lnt\k\^q^, + . . . , \k\x,qx7^0, 



Q x 

To check (formally) the validity of the expansion ( 2.3 ) we now return to the full system. Consider 



the nonlinear term e^^ in (1.2). and observe the following: 

• Case Qx everywhere on an open subinterval of [0, 2tt]. Then, on one hand, the left-hand 
side of the first equation in ( |1.2[ ) is of order at most. On the other hand, the term 
e^^ is negligible with respect to if and only if the asymptotic velocity satisfies |fc| < 1 
and is of the same order if |fc| — 1. 

• Case qx = on an open subinterval of [0,27r]. Then, e^^ is negligible with respect to 
t~'^, and no condition on the velocity k is required on that interval. 

• Case qxi^o) = at some isolated point xq. Then, no definite conclusion can be obtained 
and a "competition" between (which may approach the interval [0,1]) and qx (which 
approaches zero) is expected. 

Similarly, at least when Ifcl^-gj. ^ 0, the nonlinear term PxQx is of order \nt and, therefore, 



negligible with respect to t^d'^l (given by the left-hand side of the second equation in (1.2|) if 



and only if the asymptotic velocity is \k\ < 1. Points where \k\x or qx vanish lead to a less singular 
behavior and condition on the velocity can also be relaxed. 

The formal derivation above strongly suggests that we seek solutions to the full nonlinear equa- 



tions admitting an asymptotic expansion of the form (2.3), that is 

P^-k\nt + Lp + o(l), Q^q + t^'' i^J + o(l)), 

where fc > and Lp^q^ip are prescribed. In other words, these solutions asymptotically approach a 
solution of the VTD equations and, in consequence, such solutions will be referred to as asymptot- 
ically velocity term dominated (AVTD) solutions. 

Based on this analysis and extensive numerical experiments, it has been conjectured that asymp- 
totically as one approaches the coordinate singularity t = the function P{t,x)/lnt should approach 
some limit k = k[x), referred to as the asymptotic velocity, and that k{x) should belong to [0, 1) 
with the exception of a zero measure set of "exceptional values" . 

2.2. Gowdy equations as a second-order hyperbolic Fuchsian system. The first step in our 



(rigorous) analysis of the Gowdy equations (1.2) now is to write them as a system of second-order 
hyperbolic Fuchsian equations. These are equations of the form (written here for a scalar field in 
order to keep the notation simple) 

(2.5) D^v{t, x) + 2a{x)Dv{t, x) + b{x)v{t, x) — t^c^{t, x)d1v{t, x) + f{t, x, v, Dv, dxv), 

where w : (0, 5] x [/ — > M is the unknown defined on an interval J7 C M an interval (with S > 0). 
Here we use the symbol D := tdt, which denotes the time derivative operator and is singular at the 
origin t = 0; by D^v we mean D{Dv) = tdt{tdtv). We assume that all quantities are periodic in x 
with a periodicity domain U, and for the present applications we set U := [0, 27r]. The coefficients 
a and b are smooth and depend on the spatial coordinate x, only. The characteristic speed c has to 
satisfy certain properties at t = 0; however, since in our application we have c = 1, we do not need 
to discuss those here. All further details are explained in [3] . The left-hand side of the equation is 
called the principal part, and its right-hand side the source-term. In the course of the discussion, 
we will often abuse notation slightly and refer to the function / alone as the source-term. 
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In the first part [2, we determined leading-order behavior of solutions to (2.51 at < = 0, from 
the assumptions that it is "driven" by the principal part of the equation (in a well-defined manner) 
and, additionally, the source-term has a suitable decay property at t = 0. In particular, spatial 
derivatives are negligible, in a sense made precise in [5]. 



From the equation (2.5), we define 



uo{t,x) 



(x) i-'^(^) In t + u^^x) t-<'''> = 6, 



and 

(2.6) Ai := a + y^a^ — b, A2 :~ a — \/a 



Now, uq is regarded a prescribed (real-valued) function and, as far as our application to (1.2 1 is 
concerned, Ai and A2 are real-valued. The function uq is smooth on (0, S] x U, provided either 
a'^{x) b{x) for all x G U or else a^{x) — b{x) for all x G U. More generally, when there is a 
transition between these two regimes, the functions and w** have to be renormalized [3] in order 
to guarantee the smoothness of Uq, which we assume from now on. The function uq represents the 
leading-order terms in the ^^singular initial value problem with two-term asymptotic data and 
tij^j^ , as discussed in ^ and is referred to as a canonical two-term expansion. 



After multiplication by i^, the equations (1.2) immediately take the second-order hyperbolic 
Fuchsian form 

D^P = t'dlP + e^^'iDQf - t^e^P{d^Q)\ 

D^Q = t^dlQ - 2DPDQ + 2t^d^Pd^Q. 

The general canonical two-term expansion then reads 

P{t,x) ^P^{x)\ogt + P^^{x) + ... 

for the function P and, similarly, an expansion Qs,{x) logt -\- (5**(a;) -|- . . . for the function Q with 
prescribed data Q*, Q**. At this stage, we do not make precise statements about the (higher-order) 
remainders, yet. In any case, the theory in [3] does not apply to this system directly, due to the 
presence of the term —2DPDQ — with the exception of the cases = or = 0. Namely, 
this term does not behave as a positive power of t aX t — when we substitute P and Q by their 
canonical two-term expansions, but this is required by the theory. 



At this juncture, motivated by the formal discussion in Section 2.1 especially (2.4), we propose 
to add a term —2kDQ to the equation for Q where k is a prescribed (smooth, spatially periodic) 
function depending on x, only. This yields the system of equations 

D'P = t^dlP + e^P{DQf - ee^P{d^Q)\ 

(2.7) 

D'^Q - 2kDQ = t^dlQ - 2{k + DP)DQ + 2t^d^Pd^Q. 

Later, the function k will play the role of the asymptotic velocity mentioned before. The resulting 
system is of second-order hyperbolic Fuchsian form with two equations, corresponding to 

(2.8) A^^^ = A^^^ = 0, Af ' = 0, X^i^ ^ -2k. 



Here, the superscript determines the respective equation of the system (2.7). If we assume that k 



is a strictly positive function, as we will do in all what follows, the expected leading-order behavior 
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ai t — given by the canonical two-term expansions is 

P{t, x) = (x) log < + P„ (x) + . . . , 
^^■^^ Q{t, x) = (x) + Q„ + . . . 

One checks easily that the problem, which we had for the previous form of the equations, does 
not arise if P* — —k. Indeed, the canonical two-term expansion is consistent with the heuristics of 
the Gowdy equations above and we recover the singular initial value problem studied rigorously in 
[TBI [TO] and numerically in [T] . We only mention here without further notice that the case fc = 
with the logarithmic canonical two-term expansion for Q given by the first form of the equations 
above is covered by the following discussion. Furthermore, the case of possibly vanishing k may be 
also included via a suitable normalization of the asymptotic data [3] . 



3. Well-posedness theory for the Gowdy equations 



3.1. Reformulation of the problem. When P* = — fc, this function plays a two-fold role in (2.7|. 
On one hand, it is an asymptotic data for the function P and, on the other hand, it is a coefficient 
of the equation satisfied by the function Q. In order to keep these two roles of k separated in a first 
stage and instead of (2.7), we consider the system 

D^P = edlP + e^''{DQf - t^e'^'id^Qf, 
D^Q - 2kDQ = t^dlQ - 2(-P, + DP)DQ + 2t^d:,Pd^Q. 

Studying the singular initial value problem with two-term asymptotic data means that we search 
for solutions to (3.1 ) of the form (as t — > 0) 

P{t,x) = P^{x)\ogt + P^^x) + w^^\t,x), 
Q{t,x) = Q,(a;) + Q**(a;)t2'=(^) +iy(2)(t^a;)^ 



(3.2) 



for general asymptotic data P*, P**, Q*, Q**, and remainders w^^\ w^^-* which will be specified to 
be suitably "small" and belong to certain functional spaces; see below. After studying the well- 
posedness for this problem, we can always choose P.^, to coincide with —k and, therefore, recover 
our original Gowdy problem (2.7)-(2.9|. For simplicity in the presentation, we always assume that 
A: is a C°° function. 

In the following discussion, we write the vector-valued remainder as w := {w^^\ u'^^-'), and we fix 
some asymptotic data P*, P**, Q*, and Q**. In agreement with the notation in [3], the source-term 
operator F = (Pi, P2) is defined by 

F[w]{t,x) := f(^t,x,uo + w,D{uo + w),da:{uo + w)^. 



where uq is the vector- valued, canonical two-term expansion given by (3.2) and the asymptotic 
data. We write 

F[w]{t,x) (PiH(i,x),P2[u>](i,x)), 



and from (3.1) we obtain 



\2kf^Q,^+Dw^^^] 
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and 

+ 2 (td^P^ log t + td^P^.^. +td^w'^^^^(^{td.^Q^ +2d^kt'^^tliitQ^.^ +t'^^td^Q^^ +td^w'-'^'^ 

3.2. Properties of the source-term operator. To establish the well-posedness of the singular 
initial value problem for the Gowdy equations, we need first to derive certain decay properties of the 
source-term operator F. In [3], we introduced the spaces Xs^a,k and Xs^a,k associated with constants 
S,a > and non-negative integers k: in short, a function w belongs to Xs^a.k if each derivative 
d[.D™^w (with I + m < k) weighted by t^^(^)-°' is a bounded continuous map (0,(5] — > L'^{U). 



Here, A2 given by (2.61 is determined by the coefficients of the second-order hyperbolic Fuchsian 
equation under consideration. The corresponding norm || • \\s,a,k turns (X5q, j;,|| • Hi^Q^fe) into a 
Banach space. For any w G Xs.a,k, this norm has the form lluijl^.a^fc = supj^^Q^j Es,a,k[w\(t), where 
E5^a,k[w] ■ (0,(5] — >■ K is bounded and continuous. 

The spaces Xs^a.k and the associated maps Eg ^.k are defined analogously; only the weight of the 
/c-th spatial derivative is substituted b}|^i^^'^^^^"+^. It always follows that Xs^a,k C Xg^a.k- There 
is, in fact, no reason to assume a > to remain constant in the definition of these spaces, since 
no essentially new difficulty arises in the theory [3] when ck is a (spatially periodic) strictly positive 
function in C^{U). 

Let us introduce some further notation specific to the Gowdy equations. Let ^r^^ k ^® 



space defined as above and based on the coefficients of the first equation in (3.1) and, similarly. 



(2) 

let Xg^^ f_ be the space associated with the second equation. By definition, a vector-valued map 
w — (w'-^^w'^') belongs Xs^a,k precisely if w'^^^ G X^^l^^ j, and w'^^ G xf'^^^ j,, with a := (q;i,Q!2). 
An analogous notation is used for the spaces xf''^ ^, xf]^ j, and Xs,a,k- 



Now we are ready to state a first result about the source-term of (3.1 ). 

Lemma 3.1 (Operator F in the finite differentiability class). Fix any (5 > and any asymptotic 
data 

P,,P,,,Q,,Q,* e H'^iU) m>2. 
Suppose there exist e > and a continuous function a — (01,02) : U — > (0, 00)^ so that, at each 

X eu, 

(3.3a) ai{x) + e < mm {2{P,{x) + 2k{x)), 2{P,{x) + 1)), 
(3.3b) a2{x) + e<2{l-k{x)), 

(3.3c) ai (x) — 02 (x) > e -I- min (O, 2k{x) — l) , 

(3.3d) e < 1. 



Then, the operator F associated with the system (3.1) and the given asymptotic data maps Xg a 



into Xs^a+e,m-i and satisfies the following Lipschitz continuity condition: For each r > and for 
some constant C > ( independent of S ), 



Es,a 



+e,m — 1 



F[w]-F[w] {t)<CEs,o.^„Aw-w]{t), te(0,5] 



^The (slightly) more general definition given in [3^ involved a speed coefficient c, taken here to be identically unit. 
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for all w,w £ Br C Xg ^.m, where Br denotes the closed ball centered at the origin. 

In this lemma, since G H^{U), in particular, a standard Sobolev inequality implies that 
can be identified with a unique bounded continuous periodic function on U, and the inequality 
(3.3a) makes sense pointwise. 



Proof. Consider the expression of F given at the end of Section 3.1 Let w G Xg ^.m for some (so far 
unspecified) positive spatially dependent functions ai, a2, hence w^^^ G ^1 m ^^'^ ^ "^^1^02 m- 
By a standard Sobolev inequality (since m > 2 and the spatial dimension is 1), we get that 
F[w]{t,-) G H™'~^{U) for all t G (0,(5]. Namely, if m > 2 we can control the non-linear terms of 
F[w](t, •) in all generality for a given i > if any factor in any term of F[w]{t, •), after applying up to 
TO — 1 spatial derivatives, is an element in L°°(U) - with the exception of the TOth spatial derivative 
of w which is only required to be in L'^{U). This is guaranteed by the Sobolev inequalities. Having 
found that F[w\{t, ■) G H"'-'^iU) for all t G [0,5], it is easy to check that Fi[w] G ^1^2^+, q if 

(3.4) ai{x) + e<mm(2iP,{x) + 2k{x)),2{P4x) + l)y x e U. 

Even more, condition (3.4) implies that D^Fi[w] G Xg^^^^^ ^ for all /< to — 1. 

Considering now spatial derivatives, we have to deal with two difficulties. The first one is that 
logarithmic terms arise with each spatial derivative. We find d^D''Fi[w] G X^^^^^^ p for all I < to — 1 
and fc < TO — 2 and fc + Z < to — 1 (excluding first the case fc = to — 1, 1 = 0) provided 

(3.5) ai{x) + e < min (2{P,{x) + 2k{x)), 2{P,{x) + 1)) , 

A second difficulty arises in the case fc = to — 1, / = 0. Namely, since w G Xs^a,m (and not in 
Xs,a,m), it follows that in particular tdJ^w^"^^ f2k+a2 (^nd not ti+2fc+a2). ^ote that the function 
/3 which determines the behavior of the characteristic speeds at i = in is identically zero in 
the case of the Gowdy equations. The potentially problematic term is hence of the form AB with 

A ■.=t^'e^"e'"''\td^Q^ +2d^kt^HliitQ.,^, +t^Hd^Q.,^, +td^w''^^), 

B := ■ t'^'e^'-e'^'" {d^-^itd,Q, +2d,kt^HlntQ,, +t^Hd^Q,,) +tdl^w'^^^), 

originating from taking to — 1 spatial derivatives of To ensure d™^^Fi[w] G X^^^^^^ q, we 

need 

(3.6) ai{x) + e<{P^{x) + l) + iP4x) + 2k{x) + a2ix)), x e U. 



X G U. 



If (3.5) is satisfied, we have (for all x) 

ai{x) + e < min (2(P,(.t) + 2fc(a;)), 2(P,(a;) + 1)) 
< (P*(x) + l) + (P,(a;) + 2fc(x)) 



and, thus, (3.6) follows from (3.5 1. In conclusion, (3.5) is sufficient to guarantee that Fi[w] G 

(5,ai-f e.m — 1 ' 

Let us proceed next with the analysis of the term F2 [w] . If 

(3.7) ai{x) - a2{x) > e, a2ix) + e < 2{l - k{x)), x € U, 

(2) (2) 
then P2[w'] G ^5Q2_|_go- This inequality also implies that all time derivatives are in Xg^^_^_^g as 

before. We have to deal with the same two difficulties as before when we consider spatial derivatives 
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of F2['w]. On one hand, equality in (3.7) cannot occur due to additional logarithmic terms. On 
the other hand, we must be careful with the (m — l)-th spatial derivative of F2[w]. Here, the two 
problematic terms are of the form AB with either 

A :=a;"-i(ia,p* logt + td,p,,) + ta^w^^l 

or else 

A :=td^P^ logt + td.^P^^ + td^w^^\ 

B ■.=d^-\td^Q, + 2d.^kt^''t\ntQ^, + t^'^td^Q,^) + td^w^^\ 

The first one is under control provided ai{x) + 1 > 2k{x) + a2{x) + e, for all x G U, while for the 
second one it is sufficient to require e < 1. The claimed Lipschitz continuity condition follows from 
the above arguments. □ 

Obviously, positive functions ai and a2 and constants e > satisfying the hypothesis of 



Lemma 3.1 can exist only if k{x) < 1 for all x U (due to (3.3b|) except for special choices 



of data; cf. Lemm£3.3 below. Hence, we make the assumption that < k{x) < 1 for all x, which 



is consistent with our formal analysis in Section |2.1[ As a consistency check for the case of interest 



= — fc, let us determine under which conditions the inequalities (3.3) can be hoped to be sat 



isfied at all. For this, consider (3.3a) and (3.3c I in the "extreme" case a2 — e = 0- This leads to 



the condition < < 3/4, which shows that Lemma 3.1 does not apply within the full interval 
<k <l. 

It is interesting to note that Rendall was led to the same restriction in [T^, but its origin was 
not obvious in his approach. Here, we find that this is caused by the presence of the condition 



(3.3c) in particular which reflects the fact that w is an element of the space Xs^a,m. rather than 
of the smaller space Xs.a,m- Interestingly, we can eliminate this condition and, hence, retain the 
full interval < k < 1, when we consider the C°° -ca.se, instead of finite differentiability, as we now 
show. 

Lemma 3.2 (Operator F in the C°° class. General theory). Fix any 6 > and any asymptotic 
data 

, P^ ^ , Q ^ , Q G C {JJ^ . 

Suppose there exist a constant e > and a continuous functions a = {ai,a2) '■ U — ^ (0, cxi)^ such 
that, at each x G U, 

(3.8a) ai{x) + e < min {2{P4x) + 2k{x)), 2{P^{x) + 1)), 

(3.8b) a2(a;) + e<2(l-fc(a;)), 

(3.8c) Oii(x) — a2{x) > €. 

Then, for each integer m > 1, the operator F maps Xs^a,m into Xs,a+e,m-i cind satisfies the 
following Lipschitz continuity property: for each r > and some constant C > (independent of 
S), 

Es,a+e,rn-l[F[w] ^ F[w]]{t) < CEs,a,m[w - il>]{t), te{0,6], 

for all CXs 
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The proof is completely analogous to that of Lemma 3.1 Recall from [3, that only spaces 



without the tilde are necessary for the well-posedness theory in the C°°-case and hence that we 
obtain stronger control than in the finite differentiability case. Hence, the C°°-case does not require 
the condition (3.3c I. This has the consequence that k can have values in the whole interval (0, 1) as 
we show in detail later. In a special case, which will be of interest for the later discussion, however, 
we can relax the constraints for k even in the finite differentiability case. 

Lemma 3.3 (Operator F in the finite differentiability class. A special case). Fix any S > and 

any asymptotic data 

P,,P,,,Q,, e H'^iU), Q,^ const, m>2. 
Suppose there exist e > and a continuous function a = {0:1,02) '■ U — > (0, cxi)^ such that, at each 

X eu, 

ai{x) + e< 2{P^{x) + 2k{x)), 

a2{x) + e <2, 

ai{x) — a2{x) > e — 1, 

e < 1. 

Then, the operator F satisfies the conclusions of Lemma \3.1\ 

In the special case of constant asymptotic data = const, we can prove the required properties 
of F if the function k is any positive function in the finite differentiability case. The analogous 
result for the C°°-case can also be derived. 

3.3. Well-posedness of the singular initial value problem. Relying on Theorem 3.10 in [3], 
we now determine conditions that ensure that the singular initial value problem for the Gowdy 
equations is well-posed. Besides the properties of the source-operator F already discussed, we have 
to check the positivity of the energy dissipation matrix 

/5R(Ai-A2)+a ((3Ai)V77-77)/2 

iV:= ((3Ai)V?7-7/)/2 a td^c - d^^{Xi - X2)itc\nt) 

\ td^c~d,:^{Xi-X2){tclnt) ^{Xi - X2) + a - I - Dc/c^ 

for some well-chosen constant r/ > and for each of the two Gowdy equations. Here, we omit the 
upper indices order to simplify the notation, while 3? and 3 denote the real and imaginary part of 
a complex number. In the present paper, the characteristic speed c in p] is constant equal to 1, 
and all eigenvalues are real, so that the above matrix simplifies: 

^Ai-A2 + q; ~t]/2 

-ry/2 a -d^{Xi - X2){tlnt) 

-(9^(Ai - A2)(tlnt) Ai-Aa + a-l 





iV(i) := 


/ ai 


-77/2 


(3.9) 


[r 


ai 












In view of (2.8 1, this leads us to the matrix 




ai - 

for the first component and to the matrix 

^2A: + a2 -?y/2 
(3.10) iV(2) I -?y/2 a2 -2d^k{t\\it) 

-2dS{t\\it) 2fc + a2-l 
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for the second component. 

For the matrix N^^'> to be positive, it is necessary that ai{x) > 1 for all x € U. However, if 



—k, then condition (3.3a) in Lemma 3.1 in the finite differentiability case (or the corresponding one 
in Lemma 3.2 in the C°°-case) implies that ai{x) < 1. Hence, in the same way as in Rendall one 
does not arrive at a well-posedness result for the singular initial value problem yet. However, since 
the positivity of the energy dissipation matrix is the only part of the hypothesis in Theorem |3.10| 



of [3] which is is violated, we can use instead Theorem 3.12 of [3] to prove well-posedness of the 
"singular initial value problem with asymptotic solutions of sufficiently high order" . 

Let us quickly recapitulate the basics for this singular initial value problem which are discussed 
in detail in the first Part of this series. Consider any second-order hyperbolic Fuchsian equation of 



the form (2.5) and, for any given asymptotic data, define 

F[w] := F[w] + t^k^dliuQ + w). 

Let H be the operator which maps any source function /o = fa{t,x) (having a suitable behavior 
at t = 0) to the remainder w = v — uq, where v the unique solution of the ordinary differential 
equation 

D^v(t, x) + 2a{x)Dv{t, x) + b{x)v{t, x) = fo{t, x), 

consistent with the prescribed asymptotic data. Finally, set G := H o F . As is easily checked, w is 
the remainder of a solution of the full equations (consistent with the prescribed asymptotic data) 
if and only ii w — G[w], that is, if w is a fixed point of the map G. 
Set wi = and define the sequence 

Wj+i = G[wj], j = l,2, ... 

The convergence of this sequence to a fixed point is known for analytic data and for ordinary 
differential equations, only. Yet, the sequence (wj) has certain useful properties. 

On one hand, the residual of the second-order hyperbolic Fuchsian equation, i.e. the difference 
between the left- and the right-hand side is of higher order in t {at t — 0) if j is larger. Hence, the 
sequence satisfies the original equations at a higher order of approximation if we choose larger values 
of J. A disadvantage is that the higher we choose j, the more spatial derivatives of the asymptotic 
data we need to control. In any case, the main point for the current discussion of well-posedness 
for the Gowdy equations is the following one. 

Define the exponent 

(3.11) a a+ (j - 2)Ke, 

where a and e are the quantities introduced above and k < 1 is a constant which we can choose 
arbitrarily. If u = (P, Q) is a solution of the Gowdy equations corresponding to given asymptotic 
data, we set 

w := V — uq — Wj 

for some j > 1- Then it follows from our considerations in the first paper that the equation has a 
unique solution with remainder w G Xs,s^k for some fc, provided j is large enough so that the energy 
dissipation matrix (evaluated with a) is positive. Hence, our previous discussion implies that the 
singular initial value problem with asymptotic solutions of order j is well-posed provided one of 
the previous lemmas applies. (This is only true if the asymptotic data functions are sufficiently 
regular.) 

We can be more specific about what we mean by j being "sufficiently large" , and we now make 
some choice for the parameters ai, a2 and e, consistent with Lemma |3.2[ which will allow us to 
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estimate the required size of j. We make no particular effort to choose these quantities optimally, 
but still the goal is to choose j "reasonably" small. Henceforth, we restrict to the C°°-case and 
= —k{x) with < k{x) < 1 for all x & U. We introduce positive constants and ^2 (with 
further restrictions later) and the function xi^) '■— ^ ~ 2|a; — 1/2|. The condition (3.8a) states that 
we must choose ai{x) and e so that ai{x) + e < x(/c(a;)). We set 

(3.12) aiix) := 1 - ^ 4ik{x) - 1/2)^ + lij, 

and find xik{x)) — ai{x) > \/l + fJ^f — 1 for all x e U, provided < k{x) < 1. Similarly, we set 

(3.13) a2{x) 1 - ^4{k{x) - 1/2)^ + fij, 

and it follows that ai{x) — 02(2;) > — -^/l + /i^ for fi2 > fJ^i- For the conditions (3.8a) and 

(|3.8c|) to hold true, we have to choose 

< < H2, and < e < min ( 1 + fif - 1, J 1 + ^'^ - J 1 





Condition (3.8b) is then satisfied automatically. 

Now, assume in what follows that k{x) e (1/2 — Afc, 1/2 + Afc) for all x £ U for a constant 
Afc e (0, 1/2). Then it is clear that both functions ai and 02 are positive for all such k{x) if and 
only if 

fii<fi2< v/l-4(AA:)2. 

This assumption will be made in the following. In Theorem 3.12 of [3], we could choose j as small 
as possible if we pick the maximal allowed value for e. Hence, we set 





We find easily that 
provided 



1 < \ 1 



M? < \ilil+2^1 + f4-2), 



and check that this is consistent with the condition < ^1 < /i2 made before. In order to make a 
specific choice, we assume this inequality for fii and hence obtain that 



(3.14) 



1 



1. 



Now, in order to make the energy dissipation matrix positive, we must choose j so that for all 

X eu, 

ai{x) :— ai{x) + {j — 2)ne > 1, 

52(2;) := a2{x) + (j - 2)Ke > 1 - 2k{x); 



cf. (|3.11|). These two inequalities are satisfied for all functions k under our assumptions if in 



particular 
(3.15) 



J >2 
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(3.16) 



Ml 



m1 



2a/1 



f4 



2, 



In any case, we choose the niaximal value for fii 

2 

since this minimizes the vahie on the right side of (3.15). We find that for this value of fii, the 
right side of (3.15) is monotonically decreasing in /i2 and diverges to +00 for /Z2 — >■ for all values 
of Afc. 

Theorem 3.4 (Well-posedness theory for the Gowdy equations). Consider some asymptotic data 

= -k, Q,, g,, e C°°(C/), 

where k is a smooth function U — > (1/2 — AA;, 1/2 + Afc) for a constant Afc S (0, 1/2). Then, for 
the Gowdy equations with these prescribed data are well-posed in certain weighted Sobolev spaces. 

More precisely, the singular initial value problem with asymptotic solutions of order j has a 
unique solution with remainder w € Xs^a+{j-2)K,e,oo for some sufficiently small S > and some 



K < 1. Here, the exponents a = (ai,a2) md e are given in (3.12), (3.13), and (3.14) explicitly 
in terms of the data and parameters ^1,^12 chosen such that /ii is an explicit expression in ^2 
given in (3.16) while fi2 is a sufficiently close to (but smaller than) ^Jl — 4(Afc)^, and the order of 
differentiation j satisfies 

9 

J>2- 



3 - 4(Afc)2 + 2v^2-4(Afc)2 - 2 

The above condition implies that to reach Afc — > we need j > 7, while Afc — > 1/2 requires 
j — )■ 00. Although our estimates may not be quite optimal, the latter implication cannot be avoided. 



3.4. Fuchsian analysis for the function A. So far we have considered the equations (1.2) for 



P and Q. We can henceforth assume that these equations are solved identically for all t > (and 
t < 5 ior some 5 > 0) and that hence P and Q are given functions with leading-order behavior 
(2.9) and remainders in a given A5 „ fe. The equations which remain to be solved in order to obtain 
a solution of the full Einstein's field equations are (1.3) and (1.4). In particular we are interested 
in the function A in order to obtain the full geometrical information. We must compute A also 
as a singular initial value problem with "data" on the singularity analogously to P and Q. The 
following discussion resembles the previous one and we only discuss new aspects now. 

Clearly, the three remaining equations are overdetermined for A and hence solutions will exist 



only under certain conditions. Let us define the following "constraint quantities" from (1.4) 



Ci{t,x) := -dtK + t{P.,f^ 
C2{t,x) :=-A, 

Moreover, we define 

H{t,x) := -Afi + A^ 



2P,.DP 



+ t{dtPf + 
2e^PQ,DQ. 



2P 



from (1.3 1. From the evolution equations for P and Q, we find 
(3.17) dtCi = d,C2 + H, dtC2 = d^Ci. 



These equations have the following consequences. Suppose that we use (1.4b) as an evolution 
equation for A. This implies that Ci = for all t > 0. Moreover, suppose that we prescribe data 
at some to > (indeed to is allowed to be zero later) so that C2{to,x) = for all x € U. Then 
the equations imply that H = and C2 = for alH > and thus we have constructed a solution 
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of the full set of field equations. Alternatively, let us use (1.3) as the evolution equation for A, i.e. 
H = 0. Suppose that we prescribe data so that Ci{to,x) = C2ito,x) = at some to. It follows 
that Ci = C2 = for all i > because the evolution system (3.171 for Ci and C2 is symmetric 
hyperbolic. Again, Einstein's field equations are solved. 

Now, we want to consider the case to = 0. First note that (3.17) is regular even a,t t — 0. Suppose 
that P and Q are functions with leading-order behavior (2.9) and remainders in a given Xs,a,k with 
A: > 1. If there exists a function so that 

(3.18) K{t,x) = K{x)h\t + K^{x) +wz{t,x) 
with converging to zero in a suitable norm at i = and 

(3.19) K{x)^k^{x), A„(a;) = Ao + 2 /" fc(£)(-9sP„(£) -f 2e2^"(*)Q**(i)ajQ*(i)) di, 
where Aq is an arbitrary real constant, then 



lim C2 



0. 



Let us first use (1.4b I as an evolution equation for A. One can show easily that there exists a 
unique solution for A for t > which obeys the two-term expansion above. Our discussion before 



implies that ( 1.3 1 is solved identically for all t > 0. Hence we obtain a solution of the full Einstein's 



field equation. Alternative, choose (1.3) as the evolution equation for A now. This equation can be 



written in second-order hyperbolic Fuchsian form 

D^A - t^dlA = {td^Pf + {DA - {DPf) + e^^iitd^Qf - {DQf). 



Indeed this equation is compatible with the leading-order expansion ( |3.18 ) at i = and we can 
show well-posed of this singular initial value problem in the same way as we did for the functions 
P and Q before using the results of [3]. In particular, for any asymptotic data A* and A**, there 
exists a unique solution of this equation A with remainder in a certain space A5 q, fc. Uniqueness 



implies that the solution A of this equation coincides with the solution for A obtained using (1.4b) 
as the evolution equation. Hence, we have Ci = C2 = for all t > 0, and thus also this method 
yields a solution of the full Einstein's field equations. 

Note that periodicity implies that the asymptotic data for P and Q must satisfy the relation 



k{x){—diP^.j,{x) + 2e 



Q*^,{x)dxQ*{x)) = 



in the case of smooth solutions. 



4. Numerical solutions of the singular initial value problem 

4.1. The Fuchsian numerical algorithm. We proceed now with the numerical approximation of 
the singular initial value problem associate with second-order hyperbolic Fuchsian equations. The 



approximation algorithm proposed now is motivated by our proof of Proposition 3.4 in [3 . For linear 
source-terms, at least, we have shown that the solution of the singular initial value problem can be 
approximated by solutions to the regular initial value problem. The regular initial value problem 
is defined by data not at the singular time t — Q, but rather at some to > Oj E^nd by considering 
the evolution toward the future (i.e. away from the singular time t — Q). Then, letting to — 0, 
the sequence of solutions to these regular problems, referred to as approximate solutions, converges 
toward the solution of the singular initial value problem. In [3], we established an explicit error 
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estimate for these approximate solutions and the resuh should extend to nonlinear source-terms 
satisfying a Lipschitz continuity conditions. 

The regular initial value problem for second-order hyperbolic equations corresponds to the stan- 
dard initial value problem for a system of (non-linear) wave equations, and there exists a huge 
amount of numerical techniques for computing solutions. However, a second-order Fuchsian equa- 
tion written out with the standard time-derivative dt (instead of D) clearly involves factors \/t 
or Although these are finite for the regular initial value problem, they still can cause severe 

numerical problems when the initial time io approaches zero, due to the finite representation of 
numbers in a computer. In order to solve this problem, we introduce a new time coordinate 

r := hit, 

and observe that D ^ dr- For instance, the following Eulcr-Poisson-Darboux equation was already 
treated in [3' 

(4.1) dlv~\drv~e^^dlv^0, 

where v is the unknown and A is assumed to be a non-negative constant. Therefore, there is 
no singular term in this equation; the main price which we pay, however, is that the singularity 
t = has been "shifted to" t ~ —oo. Another disadvantage is that the characteristic speed of 
this equation (defined with respect to the r-coordinate) is e'^ and hence increases exponentially 
with time. For any explicit discretization scheme, we can thus expect that the CFL-condition will 
always be violated when we evolve into the future from some time on. We must either adapt 
the time step to the increasing characteristic speeds, or, when we decide to work with fixed time 
resolution, accept the fact that (for any given resolution) the numerical solution will eventually 
become instable. However, this is not expected to be a severe problem, since one can compute the 
numerical solution with respect to the r-variable until some finite positive time when the numerical 
solution is still stable and if necessary switch to a discretization scheme based on the original t- 
variable afterwards. All the numerical solutions presented in this paper are obtained with respect 
to the T-variable without adaption. 

We can simplify the following discussion slightly by writing (and implementing numerically) the 
equation not for the function v but for the remainder w — v — uq with uq being the canonical two- 
term expansion defined in 



to the proof of Proposition 
by initial data 



3.24) in j3j which is determined by given asymptotic data. According 



3.4 in [3], the remainder w of any approximate solution is determined 



w{to,x) = 0, drW{TQ,x) — 0, X £U, 



for some tq G M successively going to — cx). 

Inspired by Kreiss et al. in [17] and by the general idea of the "method of lines" , we proceed 
as follows to discretize the equation. First we consider second-order Fuchsian ordinary differential 
equations (written for scalar equations for simplicity) 

dlw + 2adrW + bw = /(r), 

where / is a given function and the coefficients a and b are constants. We discretize the time 
variable r so that r„ tq + nAr, w„ :— w{Tn) and /„ := /(r„) for some time step At > 0. Then 
the equation is discretized in second-order accuracy as 

(4.2) — + 2a — + 6u,„ = /„. 
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Solving this for Wn+i allows to compute the solution w at the time t„+i from the solution at the 
given and previous time t„ and t„_i, respectively. At the initial two time steps tq and ti, we set, 
consistently with the initial data for w at tq above, 

(4.3) wo = 0, w;i = i(AT)V(To). 

We will refer to this scheme as the Fuchsian ODE solver. 

The idea of the "method of lines" for Fuchsian partial differential equations is to discretize also 
the spatial domain with the spatial grid spacing Ax and to use our Fuchsian ODE solver to integrate 
one step forward in time at each spatial grid point. The source-term function /, which might now 
depend on the unknown itself and its first derivatives, is then computed from the data on the 
current or the previous time levels. Here we understand that spatial derivatives in the source-term 
are discretized by means of the second-order centered stencil using periodic boundary conditions. 
A problem is that /, besides spatial derivatives, can also involve time derivatives of the unknown 
w (in fact this can be the case for Fuchsian ordinary differential equations when the source term 
depends on the time derivative of the unknown). In order to compute those time derivatives in a 
systematic manner in second-order accuracy, i.e. without changing the stencil of the Fuchsian ODE 
solver, we made the following choice. In the code we store the numerical solution not only on two 



time levels, as it is necessary up to now for the scheme given by (4.2) and (4.3), but on a further 
third past time level. The time derivatives in the source-term can then be computed from data at 
the present and previous time steps only as follows 

drW{Tn) = ^ + 0((Ar) ). 

For this, we need to initialize three time levels and hence we set 

W2 = 2(Ar)V(To), 



in addition to (4.3 ) 



4.2. Euler-Poisson-Darboux equation. We first tested the Fuchsian ODE solver on Fuchsian 
ordinary differential equations. However, in the following, we consider the P.D.E.'s set-up directly 



and present numerical results for the Euler-Poisson-Darboux equation (4.1). The reason for con- 
sidering this equation is that can be considered as a linear version of each of the Gowdy equations. 
Recall from [3] that the singular initial value problem with two-term asymptotic data for this equa- 
tion is well-posed for < A < 2 and becomes ill-posed for A = 2. We study now this singular initial 
value problem numerically for A > 0, i.e. we look for solutions 

v{t, x) = u^,{x) + u^^,{x)t^ + w{t, x), 

with remainder w. We choose the asymptotic data — cosx, u^^ = 0. Note that in this case, this 
leading-order behavior is consistent even with the case A = 0. But according to the discussion in [3], 
it is not consistent with A = 2, and we expect that this becomes visible in the numerical solutions. 
For u** = and < A < 2, we can show that the leading-order behavior of the remainder at i = 
is 

- {-2(2^^/ + 8(2-A)(4-A) ^' + " 

First we check that the numerical solutions converge in second-order when Ar and Ax are 
changed proportionally to each other for a given choice of initial time tq > 0. We do not discuss 
this further here, but eventually we choose the resolution so that discretization errors are negligible 
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(a) A = 1.99. 
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(b) A = 1.90. 
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(c) A = 1.00. (d) A = 0.01. 

Figure 1 . Euler-Poisson-Darboux equation. 




in the following discussion; the same is true for round-off errors, as we come back to later. In 
Figure [l] we show the following results obtained with N = 20, At = 0.003. Here, N is the number 
of spatial grid points, i.e. one has Aa; = 2Tr/N. We find that for this CFL-parameter At/ Ax « 0.01, 
the runs are stable until r « 5. For each of the plots of Figure [T] we fix a value of A and study the 
convergence of the approximate solutions to the (leading-order of the) exact solution ( |4.4[ ) for various 
values of the initial time tq. We plot the value at one spatial point x — only. The convergence 
rate for tq — >■ —oo is fast if A = 1 or A = 0.01, but becomes lower, the more A approaches the value 
2, where it becomes zero. This is in exact agreement with our expectations and consistent with the 
error estimates derived in Hence the numerical results are very promising. 

To close the discussion of this test case, let us add some comment about numerical round-off 
errors. All numerical runs in this paper were done with double-precision (binary64 of IEEE 754- 
2008), where the real numbers are accurate for 16 decimal digits. However, for the case tq = 
—20 for instance, the second spatial derivative of the unknown in the equation is multiplied by 
exp(— 40) « 10~^^ at the initial time which is not resolved numerically and hence could possibly 
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T(,=-2.5, At=1.26-10", 
To=-5.0, At=1 .26-10", 
T„=-10.0, At=1.26-10, 
Tn=-10.0, At=0.63-10 



10" 
10-^ 
10-* 

lO"" 

d lo-'" 
10-^^ 
10-'* 
10-"= 



,=■2.5, At=1 .26-10-' 
,=■5.0, At=1 .26-10-, 
=■10.0, At=1 .26-10-, 
=■10.0, At=0.63-10 




(a) k = 0.5. 



(b) k = 0.9. 



Figure 2. Homogeneous pseudo-polarized case (Go-wdy equations). 

lead to a significant error. This, ho-wever, does not seem to be the case since we obtained virtually 
the same numerical solution -with quadruple precision (binary 128 of IEEE 754-2008), i.e. -when the 
numbers in the computer are represented "with 34 significant decimal digits. 

4.3. Singular initial value problem for the Go"wdy equations. We continue our discussion 
"with the singular initial value problem for the Go"wdy equations. In all of "what follo"ws we consider 
the singular initial value problem "with t"wo-term asymptotic data for the Gowdy equations. The 
fact that this "works very "well and "we get good convergence can be seen as an indication that this 
singular initial value problem is "well-posed. Recall from Theorem |3.4| that our analytical techniques 
are only sufficient to sho"w that the initial value problem "with asymptotic solutions of sufficiently 
high order is "well-posed for the Go"wdy equations. 

Test 1. Homogeneous pseudo-polarized solutions. Before we proceed "with "interesting" solutions of 
the Go"wdy problem, let us start "with a test case for "which "we can construct an explicit solution. 
Let P and Q be solutions of the polarized equations in the homogeneous case, i.e. one set Q — 
and P{t, x) = P{t)- In this case, it follo"ws directly that the exact solution of the Go"wdy equations 
is _ _ 

P{t) = -fclnt + P,*, 

"where both k and P^,^ are arbitrary constants. By a reparametrization of the Killing orbits of the 
form 

^2 = a;2/^ + X3/V2, X3 ^ -X2/V2 + X3/V2, 

"where X2, and X3 are the coordinates used to represent the orbits of the polarized solution above, 
the same solution gets reexpressed in terms of functions 

(4.5) P = lncosh(-A:lni -f P„), Q = tanh(-fclnt -I- P**). 



Of course, these functions (P, Q) are again solutions of ( 1.2 1. Asymptotically at i = 0, they satisfy 

(4.6) P = -/clnt+(P,* -ln2) + ..., Q ^ I ~ 2e-^^"t^'' + . . . , 

from "which we can read off the corresponding asymptotic data. 
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(a) A = 0.2. 



(b) A = 0.4. 



Figure 3. Convergence for general Gowdy solutions. 



Now we compute the solutions corresponding to these asymptotic data numerically and compare 
them to the exact solution (4.5|. We pick P** = 1, so that P** = 1 — ln2, Q^, = 1 and Q** = — 2e~^. 
Since the solution is spatially homogeneous - in fact this is an ODE problem ~ we only need to do the 
comparison at one spatial point. The results are presented in Figure [2] where we plot the difference 
of the numerical and the exact value of Q versus time for various values of tq. In the first plot, this 
is done for k = 0.5 and in the second plot for k = 0.9. The plots confirm nice convergence of the 
approximate solutions to the exact solution. The fact that each approximate solution diverges from 
the exact solution almost exponentially in time is a feature of the approximate solutions themselves 
and not of the numerical discretization, as is checked by comparing two different values of At in 
these plots. From our experience with the Euler-Poisson-Darboux equation, we could have expected 
that the convergence rate is lower in the case k = 0.9 than in the case k = 0.5 (note that k plays the 
same role A/2). In the case of the Euler-Poisson-Darboux equation, the rate of convergence decreases 
when A approaches 2, due to the influence of the second-spatial derivative term in the equation. In 
the spatially homogeneous case here, however, this term is zero and hence this phenomenon is not 
present. The "spikes" in Figure [2] are just a consequence of the logarithmic scale of the horizontal 
axes and the fact that the numerical and exact solutions equal for one instance of time. 

Test 2. General Gowdy equations. Now we want to study the convergence for a "generic" inhomo- 
geneous Gowdy case (still ignoring the equation for the quantity A). Here we choose the following 
asymptotic data 

k{x) = 1/2 + Acos{x), = 1.0 + sin(a;), 

P„ = 1 -ln2-Kcos(a;), Q,, = -2e^^ 

with a constant A e (—1/2, 1/2). We do not know of an explicit solution in this case. In Figure[3j we 
show the following numerical results for A = 0.2 and A = 0.4, respectively. For the given value of A, 
we compute five approximate solutions numerically with initial times tq = —30, —35, —40, —45, —50 
numerically, each with the same resolution Ar = 0.01 and N = 80. The resolution parameters have 
been chosen so that the numerical discretization errors are negligible in the plots of Figure |3] Then, 
for each time step for r > —30, we compute the supremum norm in space of the difference of the 
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remainders w^^^ of the two approximate solutions given by tq — —30 and tq = —35. In this way we 
obtain the first curve in each of the plots of Figure [3j The same is done for the difference between 
the cases tq = —35 and tq = —40 for all t > —35 to obtain the second curve etc. Hence these 
curves yield a measure of the convergence rate of the approximation scheme, without referring to 
the exact solution. In agreement with our observation for the Euler-Poisson-Darboux equation, the 
convergence rate is high if k is close to 1/2 and becomes lower, the more k touches the "extreme" 
values fc = and k = 1. 

Much in the same way as for the Euler-Poisson-Darboux equation we find that double precision 
is sufficient for these computations despite of the fact that exp(2r) is 10"'*'* for r = —50. 

4.4. Gowdy spacetimes containing a Cauchy horizon. The papers [HI HOI [HI ttH H] were 
devoted to the construction and characterization of Gowdy solutions with Cauchy horizons, in 
particular in order to prove the strong cosmic censorship conjecture in this class of spacetimes. 
Spacetimes with Cauchy horizons are expected to have saddle and physically "undesired" properties, 
in particular they often allow various inequivalent smooth extensions, i.e. the Cauchy problem of 
Einstein's field equations does not select one of them uniquely. Some explicit examples are known, 
but most of the analysis is on the level of existence proofs and asymptotic expansions. 

Hence, it is of interest to construct such solutions numerically and analyze them in much greater 
detail than possible with purely analytic methods. Constructing these solutions numerically, how- 
ever, is delicate since the strong cosmic censorship conjecture suggests that they are instable under 
generic perturbations. It can hence often be expected that numerical errors would most likely "de- 
stroy the Cauchy horizon" . This is so, in particular, when the singular time at t = is approached 
backwards in time from some regular Cauchy surface at i > 0. 

In the Gowdy case, however, there are clear criteria for the asymptotic data so that the corre- 
sponding solution of the singular initial value problem has a Cauchy horizon (or only pieces thereof; 
cf. below) at t = 0, as discussed in '9' for the polarized case and in QAj for the general case. Our 
novel method here allows us to construct such solutions with arbitrary accuracy and it can hence be 
expected that this allows us to study the saddle properties of such solutions. Our main aim for the 
following is to compute such a solution and with this demonstrate the feasibility of our approach. 
A follow-up work will be devoted to the numerical construction and detailed analysis of relevant 
classes of such solutions. 

Motivated by the results in [3], we choose the asymptotic data as follows 

J 1, X € [7r,27r], 

g*(x) = 0, 

A*(x) = k'^ix), 

With these asymptotic data, the corresponding solution has a smooth Cauchy horizon at {t, x) G 
{0} X (tt, 27r) (namely where k = 1), and a curvature singularity at {t,x) € {0} x (0,7r) (namely 
where < fc < 1). Note that the function k is smooth everywhere (but not analytic). Our analysis 
in Section [3] shows that we are allowed to set fc = 1 at some points since dxQ* = 0. This motivates 
our choice of Q^. With this, our choice of Q** implies that the solution is polarized on the "domain 



* ( X ) 



1/2, 



A„(a;) = 2 



0, x€ [tt, 27r], 

^-i/x^-i/iTT-x)^ a;e(0,7r). 
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Figure 4. Convergence analysis for spacetimes with a Cauchy horizon. 

of dependence"!^ of the "initial data" interval (tt, 2tt). All data were chosen as simple as possible to 
be consistent with the constraints. 

First we repeated the same error analysis as for the previous Gowdy case, see Figure |4] For all 
the runs in the plots, we choose N — 500, At = 0.005 which guarantees that discretization errors 
are negligible in the plot. We find that our numerical method allows us to compute the Gowdy 
solution very accurately. Here, we solve the full system for {P,Q,A). 

In Figure[5) we show the numerical solution obtained from N — 1000, At = 0.0025 and To — —18. 
We plot the Kretschmann scalar at two times t = —10 and t = 0. Hence, near the time t = 
(corresponding to t = — oo), the Kretschmann scalar is large on the spatial interval (0,7r) while 
it stays bounded at (tt, 27r). At the later time, the curvature becomes smaller as expected. We 
also plot the remainders w^^^ and w^^^ of P and Q, respectively. It is instructive to study how the 
polarized region inside (tt, 27r) gets "displaced" by the non-polarized solution. We refer the reader 
to [1] for further investigations, especially a study of past-directed causal geodesies approaching the 
t = 0-hypersurface near the boundary point x = n at the intersection of the Cauchy horizon and 
the curvature singularity. Furthermore, in [4] we will discuss trapped surfaces in a neighborhood of 
t = 0. 

Acknowledgements 

The authors were partially supported by the Agence Nationale de la Recherche (ANR) through 
the grant 06-2-134423 entitled "Mathematical Methods in General Relativity" (MATH-GR). A first 
draft of this paper was written during the year 2008-2009 when the first author (F.B.) was an ANR 
postdoctoral fellow at the Laboratoire J.-L. Lions. 

References 

[1] Amorim P., Bernardi C, and LeFloch P.G., Computing Gowdy spacetimes via spectral evolution in future and 

past directions, Class. Quantum Grav. 26 (2009), 1-18. 
[2] Andersson L. and Rendall A.D., Quiescent cosmological singularities, Commun. Math. Phys. 218 (2001), 479-511. 
[3] Beyer F. and LeFloch P.G., Second-order hyperbolic Fuchsian systems. General theory, arXiv: 1004.4885. 



The notion of "domain of dependence" for the singular initial value problem follows easily from the energy 
estimate established in [3]. 



SECOND-ORDER HYPERBOLIC FUCHSIAN SYSTEMS. 



23 




Figure 5. Behavior of solutions with a Cauchy horizon. 



[4] Beyer F. and LeFloch P.G., in preparation. 

[5] Choquet-Bruhat Y. and Isenberg J., Half-polarized U(l) symmetric vacuum spacetimes with AVTD behavior, J. 

Geom. Phys. 56 (2006), 1199-1214. 
[6] Choquet-Bruhat Y., Puchsian partial differential equations, "WASCOM 2007" — 14th Conference on Waves and 

Stability in Continuous Media, World Sci. Publ., Hackensack, NJ, 2008, pp. 153-161. 
[7] Choquet-Bruhat Y., General relativity and the Einstein equations, Oxford Math. Monographs, Oxford Univ. 

Press, Oxford, 2009. 

[8] Chrusciel P., On spacetimes with U{1) x U{1) symmetric compact Cauchy surfa<;es, Ann. Phys. 202 (1990), 
100-150. 

[9] Chrusciel P., Isenberg J., and Moncrief V., Strong cosmic censorship in polarized Gowdy spacetimes. Class. 
Quantum Grav. 7 (1990), 1671-1680. 

[10] Chrusciel P. and Isenberg J., Nonisometric vacuum extensions of vacuum maximal globally hyperbolic space- 
times, Phys. Rev. D 48 (1993), 1616-1628. 

[11] Chrusciel P. and Lake K., Cauchy horizons in Gowdy spacestimes. Class. Quantum Grav. 21 (2004), S153-S169. 

[12] Eardley D. and Moncrief V., The global existence problem and cosmic censorship in general relativity. Gen. 
Relat. Grav. 13 (1981), 887-892. 



24 FLORIAN BEYER AND PHILIPPE G. LbFLOCH 

[13] Gowdy R.H., Vacuum space-times with two parameter spacelike isometry groups and compact invariant hyper- 
surfaces: Topologies and boundary conditions, Ann. Phys. 83 (1974), 203—241. 

[14] Hennig J. and Ansorg M., Regularity of Cauchy horizons in S'^ X Gowdy spacetimes. Class. Quantum Grav. 
27 (2010), 065010. 

[15] Iscnbcrg ,J. and Moncricf V., Asymptotic behavior of the gravitational field and the nature of singularities in 

Gowdy spacetimes, Ann. Phys. 99 (1990), 84-122. 
[16] Kiclienassamy S. and Rondall A.D., Analytic description of singularities in Gowdy spacetimes. Class. Quantum 

Grav. 15 (1998), 1339-1355. 

[17] Kreiss H.O., Petersson N.A., and Ystrom J., Difference approximations for the second-order wave equation, 

Siam J. Numer. Anal 40 (2002), 1940-1967. 
[18] Moncrief v., Global properties of Gowdy spacetimes with X R topology, Ann. Phys. 132 (1981), 87-107. 
[19] Rendall A.D., Puchsian analysis of singularities in Gowdy spacetimes beyond analyticity. Class. Quantum Grav. 

17 (2000), 3305-3316. 

[20] Ringstrom H., Asymptotic expansions close to the singularity in Gowdy spacetimes. Class. Quantum Grav. 21 

(2004) , S305-S322. 

[21] Ringstrom H., Curvature blow up on a dense subset of the singularity in T3-Gowdy, J. Hyperbolic DijJ. Eqs. 2 

(2005) , 547-564. 

[22] Ringstrom H., Strong cosmic censorship in T^-Gowdy spacetimes, Ann. of Math. 170 (2009), 1181-1240. 



